In [370]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from skimage import io
from skimage.viewer import ImageViewer
from numpy import linalg
from matplotlib import cm as cm
In [371]:
egImage = "https://upload.wikimedia.org/wikipedia/commons/7/7d/Dog_face.png"
img = io.imread(egImage, as_grey = True)
plt.imshow(img, cmap = cm.gray)
plt.show()
In [372]:
print (type(img))
print (str(img.shape))
In [373]:
img2 = img[5:70,0:70]
plt.imshow(img2, cmap = cm.gray)
plt.show()
In [374]:
?imread
In [380]:
img.shape
Out[380]:
In [381]:
img2.shape
Out[381]:
In [382]:
?linalg.svd #Full SVD reconstruction means diagonlization? Reduced SVD reconstruction means SVD?
In [384]:
u, s, vt = linalg.svd(img, full_matrices=True) #img es simetrica
print(u.shape)
S = np.diag(s)
print(S.shape)
print(vt.shape)
In [385]:
S=np.diag(s)
img_rec = np.matmul(np.matmul(u,S),vt)
print(img_rec.shape)
In [386]:
u, s, vt = linalg.svd(img2, full_matrices=False) #img2 NO es simetrica
print(u.shape)
S = np.diag(s)
print(S.shape)
print(vt.shape)
In [387]:
img2_rec = np.matmul(np.matmul(u,S),vt)
print(img2_rec.shape)
Y, graficamos la aproximación SVD...
In [388]:
plt.imshow(img2_rec, cmap = cm.gray)
plt.show()
In [389]:
k=1000
In [390]:
u_k = u[:,:k]
S_k = S[:k,:k]
vt_k = vt[:k,:]
img_k = np.matmul(u_k,np.matmul(S_k,vt_k))
plt.imshow(img_k, cmap = cm.gray)
plt.show()
La SVD sirve para almacenar de forma mas compacta imagenes ya que la matriz de aproximación de grado k contiene la información necesaria para replicar la imagen original
In [398]:
#Primero hagamosla para una matriz particular y despues veamos el tema de la generalización. Sea A una matriz cualquiera:
A = np.array([[1,2,3],[4,5,6],[7,8,9]])
Au, As, Avt = linalg.svd(A, full_matrices=True)
#As = np.diag(As)
A_inv = np.linalg.inv(A)
Au_inv = np.linalg.inv(Au)
As_inv = np.diag(np.array([1/x if x != 0 else 0 for x in As]))
Avt_inv = np.linalg.inv(Avt)
print(' A_inv')
print(A_inv)
print(' Au_inv')
print(Au_inv)
print(' As_inv')
print(As_inv)
print(' Avt_inv')
print(Avt_inv)
In [399]:
#Verificacion del calculo de las inversas
print(' Au y Au inversa')
print(np.dot(Au,Au_inv))
print(' As y As inversa')
print(np.dot(As,As_inv))
print(' Avt y Avt inversa')
print(np.dot(Avt,Avt_inv))
La definición de la psuedo inversa es la siguiente. Partiendo de la SVD $ \mathbf A = U \Sigma V^{T} $ por lo que la pseudoinversa de A será $ A^{+}=V \Sigma^{-1} U^{T} $
In [400]:
#NO ME DA PREGUNTAR MAÑANA
A_pseudo_inv = np.dot(Avt_inv,np.dot(As_inv,Au_inv))
print(" Metodo propio")
print (A_pseudo_inv)
print(" Metodo numpy ")
print(A_inv)
La función:
In [395]:
def get_pseudoinverse(matriz):
"Regresa la pseudo inversa de la matriz dada"
Au, As, Avt = linalg.svd(matriz, full_matrices=True)
#As = np.diag(As)
Au_inv = Au.transpose()
As_inv = np.diag(np.array([1/x if x!=0 else 0 for x in As]))
Avt_inv = Avt.transpose()
pseudo_matriz = np.dot(Avt_inv,np.dot(As_inv, Au_inv))
return pseudo_matriz
Verifiquemos que sirve...
In [396]:
#Verifiquemos que sirve
A = np.array([[1,2,4],[1,4,6],[2,1,6]])
Ainv = get_pseudoinverse(A)
iden = np.dot(A,Ainv)
print("La matriz A por su pseudoinversa debe ser la identidad")
iden
Out[396]:
Función para solucionar sistemas de ecuaciones lineales de la forma: $\mathbf b=A\overrightarrow{x} $ usando la pseudoinversa de A talque $ \mathbf A^{+}b=\overrightarrow{x}$
In [365]:
def resuelve_ec(A,b):
#Primero obtengamos la pseudoinversa con el metodo get_pseudoinverse recien creado
A_inv = get_pseudoinverse(A)
x = np.dot(A_inv, b)
return x
Verifiquemos que sirve...
In [366]:
A = np.array([[1,2],[2,1]])
b = np.array([3,3])
print('\n matriz A =\n',A,'\n b=',b)
x = resuelve_ec(A,b)
print('\n x = ',x)
y = np.dot(A,x)
print('\n Ax = ',y)
Jugar con el sistema $ \mathbf Ax=b $ donde $ \mathbf A = [[1,1],[0,0]] $ y $ \mathbf b $ puede tomar distintos valores. ¿Cúal es la imagen de A?¿Qué pasa con la solución a la ecuación si $ \mathbf b $ no esta en la imagen de A, la solución resultante es única?
Por lo que la solución al sistema consiste en encontrar los valores de $\mathbf x_{1}, x_{2} $ para los cuales se cumple que:
</p> $$ \begin{pmatrix} a, b \\ c, d \end{pmatrix} \begin{pmatrix} x_{1} \\ x_{2} \end{pmatrix} = \begin{pmatrix} b_{1} \\ b_{2} \end{pmatrix} $$ $$ \rightarrow ax_{1}+bx_{2} = b_{1} \cap cx_{1}+dx_{2} = b_{2} $$
Por lo que si $ \mathbf c=0,d=0,b_{2}=1 $ $ \mathbf b \notin Img(A) $
Pues tendría suceder:
$$ \mathbf x_{1}+x_{2}=1 $$y
$$ 0=1 $$ $$ \rightarrow \mathbf Img(A) =\left \{ {y \in \mathbb{R}^{2} : y_{1} \in \mathbb{R} \cap 1 = 0} \right \} $$Si $\mathbf b $ no se encuentra en la $\mathbf Img(A) $ es necesario encontrar una proyección de este vector en el espacio generado por las columnas de la matriz $\mathbf A$, generalmente se usa a la matriz de proyección:
$$\mathbf P_{x} = x \left[ xx^{T} \right]^{-1} x^{T} $$Lo que implica que:
$$ \widehat{b} \mathbf = AP_{x} = A \left( x \left[ xx^{T} \right]^{-1} x^{T} \right) $$Para resolver este sistema se utiliza el método de mínimos cuadrados y la solución del método es única
Matriz A con determinante igual a 0:
En este caso la solución es {$ \mathbf x \in \mathbb{R}^{2} : x_{1}+x_{2}=1 \cap $ ¡$ 1=0 $! $ $}
Que significa que el vector $ \mathbf b \notin Img(A) $
In [367]:
# Usemos nuestras funciones para jugar con distintos valores de b
A = np.array([[1,1],[0,0]])
b = np.array([1,1])
x = resuelve_ec(A,b)
print(x)
Matriz A con determinante distinto de 0:
En este caso la solución es {$ \mathbf x \in \mathbb{R}^{2} : x_{1}+x_{2}=1 \cap 1=e^{-32}x_{2} $}
Que significa que el vector $ \mathbf b \in Img(A) $
In [368]:
A = np.array([[1,1],[0,1e-32]])
b = np.array([1,1])
y = resuelve_ec(A,b)
x = np.linalg.solve(A,b)
print(y)
print(x)
In [263]:
import pandas as pd
from sklearn.linear_model import LinearRegression
base = pd.read_csv("study_vs_sat.csv", header=0)
print(base.head())
print(base.dtypes)
Se comienza por postular a la forma del valor esperado del sat_score $ \mathbf y $ condicional a las study_hours $ \mathbf x $ como una relación lineal (MPL) en el cual si se puede asumir que $ \mathbf E(\vec{u} \mid \vec{x})=0 $ , se tiene:
$$ \vec{y} = \vec{ \beta } \vec{x} + \vec{u} \Rightarrow E( \vec{y} \mid \vec{x} ) = \vec{ \beta} \vec{x} $$Con:
$$ \vec{ \beta } = (1, \beta_{1}, ... , \beta_{n}) $$Por lo que, el problema que se desea resolver es el siguiente:
$$ \mathbf min_{ \vec{ \beta }} E(UU_{T}) = min_{ \vec{ \beta }} E(( \vec{y} - E( \vec{y} \mid \vec{x}))^2) $$Con función objetivo: $ \mathbf E(( \vec{y}- \vec{ \beta } \vec{x})( \vec{y}- \vec{ \beta } \vec{x})^{T}) $
El gradiente de este problema se encuentra al derivar con respecto a $ \vec{ \beta } $ la funcion objetivo:
$$ \bigtriangledown_{ \vec{ \beta } } = \frac{\partial E(( \vec{y}- \vec{ \beta } \vec{x})( \vec{y}- \vec{ \beta } \vec{x})^{T}) }{\partial \vec{ \beta }} $$La C.P.O de este problema obligan a que en el optimo se cumpla :
$$ \bigtriangledown_{ \vec{ \beta } } = \vec{0} $$La libreria sklearn ya contiene el modulo de regresion lineal que permite encontrar los coeficientes buscados a travéz del método fit
In [271]:
base_reg = np.matrix(base)
X, Y = base_reg[:,0], base_reg[:,1]
mdl = LinearRegression().fit(X, Y)
m = mdl.coef_[0]
b = mdl.intercept_
print("E(Y|X): E(y|x) = {0}E(x) + {1}".format(m,b))
Construye una función que reciba valores de alpha, beta y horas de estudio para devolver en vector de numpy con la predección del modelo de mínimos cuadrados correspondiente
In [273]:
def MejorPredictorLineal(alpha, beta, x):
"Regresa el resultado de la ecuacion postulada: Mejor Predictor Lineal"
X = np.array(x)
y_arr = b + m*X
return y_arr
Veamos que sucede si hay obtenemos una muestra distinta de individuos y les aplicamos el Mejor Predictor Lineal que postulamos:
In [278]:
ListaDeNuevasCalificaciones = np.array([[6, 7, 8, 9, 7]])
ResultadosPredecidos = MejorPredictorLineal(353.16, 25.33, ListaDeNuevasCalificaciones)
print('\n Resultados Predecidos E(Y|X) =', ResultadosPredecidos)
In [314]:
import itertools
import operator
print(base_reg[:,0].size)
horas_est = np.array([base_reg[:5,0]])
print(horas_est)
item = print(horas_est.size)
aux = horas_est.tolist()
merged = list(itertools.chain(*aux))
print(merged)
b = np.array([[1, 1, 1, 1, 1],aux])
print(b)
b
Out[314]:
In [315]:
sum(aux,[])
Out[315]:
In [316]:
base_reg_1 = np.matrix(base)+ b +1
In [319]:
alpha, beta = base_reg[:,0], base_reg[:,1]
In [ ]: